SIMULATIONS OF CAVITATING CRYOGENIC 

INDUCERS 


Ashvin Hosangadi/ Vineet Ahuja* and Ronald J. Ungewitter* 
Combustion Research and Flow Technology, Inc. (CRAFT Tech) 
6210 Keller’s Church Road, Pipersville, PA 18947 
Phone: 215-766-1520/ Fax: 215-766-1524 


Simulations of cavitating turbopump inducers at their design flow rate are presented. 
Results over a broad range of N„ numbers extending from single-phase flow conditions 
through the critical head break down point are discussed. The flow characteristics and 
performance of a subscale geometry designed for water testing are compared with the full- 
scale configuration that employs LOX. In particular, thermal depression effects arising 
from cavitation in cryogenic fluids are identified and their impact on the suction 
performance of the inducer quantified. The simulations have been performed using the 
CRUNCH CFD® code that has a generalized multi-element unstructured framework suitable 
for turbomachinery applications. An advanced multi-phase formulation for cryogenic fluids 
that models temperature depression and real fluid property variations is employed. The 
formulation has been extensively validated for both liquid nitrogen and liquid hydrogen by 
simulating the experiments of Hord on hydrofoils; excellent estimates of the leading edge 
temperature and pressure depression were obtained while the comparisons in the cavity 
closure region wore reasonable. 


I. Introduction 

Our focus is on the development of a computational framework to simulate cavitating liquid rocket 
turbomachinery that employ cryogenic working fluids. Liquid rocket systems are a subset of a broader class of 
pumps (e.g. refrigerant systems, boiler feed pumps, etc) where the operating temperature is elevated relative to the 
critical temperature of the fluid and thermodynamic effects of cavitation play an important role. At these operating 
temperatures, the ratio of liquid to vapor density is lower and consequently more liquid mass has to vaporize to 
sustain a cavity. Therefore evaporative cooling effects are more pronounced and result in the lowering of the mean 
fluid temperature in the cavitating region. Since the fluid thermodynamics properties (i.e. vapor pressure, density) 
are a strong function of temperature at these conditions, thermal effects suppress cavitation and lower the cavity 
pressure in a mean sense. Typically this results in improved mean performance of cryogenic pumps; liquid 
hydrogen systems being an extreme example where the pump may continue to generate head even when the fluid is 
boiling at the inlet. 

The thermal effects of cavitation were studied extensively by numerous researchers through the 1970’ s 
including: Stahl and Stepanoff 2 , Ruggeri and Moore 3 , Hord 1 , Holl 4 , and Brennen 5 among others. Stahl and 
Stepanoff 2 were the first to estimate head depression (A H v ) values due to thermodynamic effects using the so-called 
‘B -factor’ method based on a quasi-static theory where the temperature depression was estimated in terms of the 
ratio of the vapor volume to liquid volume. They provided a graph to evaluate NPSH corrections for hydrocarbons 
based on this methodology. More elaborate correlations, which included dynamic effects were given later by 
Ruggeri and Moore 3 , Hord 1 , and Holl 4 . They collected extensive experimental data of cavity pressure and 
temperature depressions for a variety of model shapes and fluids and correlated the results using variants of the B- 
factor theory. The semi-empirical procedures outlined by Ruggeri and Moore 3 continues to be used as an 
engineering tool for predicting the thermodynamic depression in pumps. A more rigorous numerical procedure was 
developed by Cooper 6 where a baratropic equation of state was used to define the two-phase mixture and thermal 
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effects were evaluated with a resulting non-dimensional vaporization parameter. Most of these techniques, however, 
require some degree of empiricism. Therefore, from a more fundamental modeling perspective, this discussion 
highlights the need for a generalized compressible formulation that takes the energy balance into account when 
simulating cavitation for cryogenic flows. 

For turbomachinery simulations, CFD technology is currently limited to simulating mean cavitating performance 
at design conditions for idealized liquids (no thermodynamic effects). Typical simulations show comparisons with 
data at design conditions for the head coefficient and the critical number at which performance breakdown 
occurs (Hosangadi 7 , Athavale and Singhal 8 , Dupont and Okamura 9 , Medvitz 10 ). However flows at off-design 
conditions, where large scale unsteadiness and high dynamic pressure loads are observed, cannot at this point be 
reliably predicted Simulation of cavitation instabilities and rotational cavitation modes in pumps have not been 
simulated by any group to the best of our knowledge. Thus, to develop a CFD simulation framework that can 
eventually simulate unsteady, off-design performance of liquid rocket turbopumps that use cryogenic working 
fluids, significant development will be required from the current level of technology. 

The focus of this paper is on the simulation of cavitating, cryogenic inducers where real fluid property variations 
have a significant impact on the suction performance of turbopumps. To simulate this class of flows, a generalized 
multi-phase formulation has been developed that rigorously models thermal effects of phase change and the 
accompanying property variations. Thermal equilibrium is assumed and fluid thermodynamic properties are 
specified along the saturation line using the NIST-12 databank (Hosangadi and Ahuja 11 ). This framework will be 
the foundation for the eventual goal of our effort to model unsteady performance of liquid rocket pumps. With this 
is mind, the equations are cast in an acoustically accurate form and is an extension of earlier work by the authors for 
idealized liquids (Ahuja 12 ). The rigorous acoustic treatment of the gas/liquid interface captures the coupling of the 
phase change process with the dynamic pressure field which is an important attribute for unsteady cloud cavitation 
problems where high amplitude, localized pressure spikes can occur (Hosangadi I3 ). 

The numerical code utilized in our simulations is the CRUNCH CFD® code, which is a multi-element based 
unstructured code 14,15 . The underlying philosophy in the CRUNCH CFD® code is to tailor the grid topology to 
resolving the dominant flow phenomena and the structural complexity of the problem. This is achieved by utilizing a 
combination of hexahedraf tetrahedral, prismatic and pyramidal elements in mesh construction. Such a framework 
is particularly attractive for complex turbomachine configurations, since high quality grids can be generated very 
efficiently with minimum skewness. 

The computational framework described above has been applied to a full-scale main LOX pump inducer. 
Furthermore, to identify the impact of thermal effects on the inducer suction performance, a sub-scale configuration 
designed for water testing has also been simulated and compared with LOX configuration. The performance for 
both configurations is simulated at the design flow rate over a range of inlet pressures ranging from essentially 
single-phase flow all the way through head breakdown at the critical N^ number. Detailed analysis of the flowfield 
over this range of N ss numbers is undertaken; this includes ascertaining the location and shape of the vapor cloud, 
the detailed pressure distributions on the blade, and the velocity distributions both upstream and downstream of the 
inducer. In particular, the additional flow losses resulting under cavitating conditions and their effect on the velocity 
profile is identified. The effect of thermal depression on the breakdown N M number is identified. We note that since 
all calculations presented here have been computed at the design flow rate large scale flow separation upstream of 
the inlet leading edge is not present and the total pressure rise is relatively steady. 


II. Multi-phase Equation System 

The multiphase equation system is written in vector form as: 

dQ dE dF dG 

— + + + = 5 + D 

3t dx dy dz 


( 1 ) 


Here Q is the vector of dependent variables, E, F and G are the flux vectors, S the source terms and D v represents 
the viscous fluxes. The viscous fluxes are given by the standard full compressible form of Navier Stokes equations. 
The vectors Q, E and S are given below with a detailed discussion on the details of the cavitation source terms to 
follow later: 
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Here, p m and h m are the mixture density and enthalpy respectively, and <j> g is the volume fraction or porosity of 
the vapor phase. The mixture energy equation has been formulated with the assumption that the contribution of the 
pressure work on the mixture energy is negligible which is a reasonable assumption for this flow regime. The 
source term for the vapor phase arises from rate of vapor mass generation due to cavitation m, and the 
corresponding source term for the energy equation ^ given as m l *h fg where h fg is the change in enthalpy resulting 

from the phase change and is a function of the local fluid temperature. 

The mixture density and gas porosity are related by the following relations locally in a given cell volume: 


Pm ~ Pgfig + PL&L ( 3 ) 

1 = ^ + (4) 

where p g , are the physical material densities of the gas and liquid phase respectively and in general are functions 
of both the local temperature and pressure. 

Thus far we have not made any statements defining the temperatures characterizing the liquid and vapor. In 
general, the liquid and vapor may not be in equilibrium locally and can have independent temperatures. Examination 
of temperature and pressure data for cavitation in Freon by Ruggeri 16 reveals that the saturation vapor pressure 
corresponding to the local fluid temperature in fact matches the local pressure measurement. This indicates local 
thermodynamic equilibrium that is exploited to significantly simplify Eqn. (1). The thermodynamic properties of 
the liquid and vapor in the cavity may now be defined by a single variable; the saturation temperature T yai . Hence, 
all thermodynamic properties (density, vapor pressure, viscosity, etc) of both the liquid and the vapor phase may be 
generated as a tabular function of the saturation temperature. In our study here, these properties were generated 
from the Standard thermodynamic database 12 available from NIST for pure fluids. 

The equation system as formulated in Eqn. (1) is very stiff since the variations in density are much smaller than 
the corresponding changes in pressure. Therefore to devise an efficient numerical procedure we wish to transform 
Eqn. (1) to a pressure based form where pressure rather than density is the variable solved for. An acoustically 
accurate two-phase form of Eqn. (1) is first derived, followed by a second step of time-scaling or preconditioning to 
obtain a well-conditioned system. We begin by defining the acoustic form of density differential for the individual 
gas and liquid phase as follows: 


1 


1 


dp g =—dP, dpi = - dP 

Cg Cl 


(5) 


Here c g is the isothermal speed of sound 


( cP ^ 

w. 


in the pure gas phase, and c L is the corresponding isothermal 


speed of sound in the liquid phase, which is a finite- value. We note that in Eqn. (5) the variation of the density with 
temperature has been neglected in the differential form. This assumption was motivated by the fact that the 
temperature changes are primarily due to the source term and not by the pressure work on the fluid i.e. the energy 
equation is a scalar equation. This simplifies the matrix algebra for the upwind flux formulation significantly, at the 
potential expense of numerical stability in a time-marching procedure. However, more importantly, there is no 
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impact on the accuracy since the fluid properties themselves are taken directly from the thermodynamic data bank 
for each fluid. 

Following the discussion above, the differential form of the mixture density p m using Eqn. (5) is written as. 


dp m ={p t ~ Pi) d</> t +—jdP 

C 4 



Here, c* is a variable defined for convenience and is not the acoustic speed. Cm, in the mixture, which will be defined 
later. Using Eqn. (6), Eqn. (1) may be rewritten as: 
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( 7 ) 


and, 

Q v = [p,u,v,w,<fig,k,£f (8) 

The numerical characteristics of the Eqn. (7) are studied by obtaining the eigenvalues of the matrix. 

The eigenvalues of the system are derived to be: 

A=(u+c m , u-c m> u,u,u,u,u) (9) 

where c m turns out to be the well-known, harmonic expression for the speed of sound in a two-phase mixture and is 
given as: 



1 

2 Pm 
c m 

The behavior of the two-phase speed of sound indicates that at either limit the pure single-phase acoustic speed 
is recovered. However, away from the single-phase limits, the acoustic speed rapidly drops below either limit value 
and remains at the low-level in most of the mixture regime. As a consequence, the local Mach number in the 
interface region can be large even in low speed flows. 

To obtain an efficient time-marching numerical scheme, preconditioning is now applied to the system in Eqn. 
(7), in order to rescale the eigenvalues of the system so that the acoustic speeds are of the same order of magnitude 
as the local convective velocities 


2 2 
Pifig Pl c l 


( 10 ) 


III. Cavitation Source Terms 

In the present effort, the cavitation source term is defined via a simplified non-equilibrium, finite rate form as 
follows: 


m , = K /PJl + K b p t <f> t 


(ID 


where the constant Kf is the rate constant for vapor being generated from liquid in a region where the local pressure 
is less than the vapor pressure. Conversely, Kb is the rate constant for reconversion of vapor back to liquid in 
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regions where the pressure exceeds the vapor pressure. Here, the rate constants are specified using the form given 
by Merkle 17 . 
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( 12 ) 


We note that for steady attached cavitation this simplified form may be adequate since the cavitation time scales 
do not interact with the fluid time scales if the cavitation rate constants are fast enough. For unsteady cavitation 
modeling, however, it becomes essential to integrate bubble dynamics within a dense cloud framework wherein both 
the number density and mean local radius of the bubbles in an evolving cloud are tracked. The development of a 
more rigorous non-equilibrium source term model is a topic of ongoing research. 

IV. Validation Study For Cavitation In Cry ogenic Fluids 

The multi-phase formulation described above has been validated extensively by simulating experiments by 
Hord 1 on a cavitating hydrofoil. Hord performed sub-scale tests using both liquid nitrogen and hydrogen in a blow- 
down tunnel. The details of the tunnel and the hydrofoil geometry are given in Fig. 1. The tunnel width is 1 inch 



Tunnel Geometry 


Quarter Caliber Hydrofoil with a Tapered Diffuser 


Figure 1. Details of Experimental Set-up for Cavitating Hydrofoil (Taken from Hord, 1973). 


while the hydrofoil width is 0.312 inches. Hence considerable 
blockage effects from the tunnel w>all are present and this necessitated 
modeling the tunnel geometry in the simulations. To ensure that the 
tunnel blockage effects were being correctly modeled we performed 
single-phase, non-cavitating simulations and have compared it to 
Hord’s non-cavitating data in Fig. 2. Excellent agreement is obtained 
giving confidence that the tunnel interaction is being captured. We 
note that the single-phase solution is insensitive to the Reynolds 
number since both liquid hydrogen and nitrogen show identical 
pressure profiles. 

Cavitating results in liquid nitrogen for one flow condition are 
discussed here. We refer the reader to Hosangadi and Ahuja 11 for 
additional validation studies in both liquid nitrogen and hydrogen. The 
operating range of liquid nitrogen varies from roughly 70K -100K. 
The variation of properties along the saturation line for liquid and 
vapor densities as well as the vapor pressure are shown in Fig. 3. 

At 89 K, the slope of vapor pressure curve indicates a 16 KPa 
increase for a 1 K change in temperature. 

The general characteristics of the cavitating flowTield are shown 
in Figure 4(a-d) for tunnel conditions in Run 290C (Table I). 



Figure 2. Non-Cavitating Pressure Distribution 
on Hydrofoil With Tunnel Blockage Modeled 
(Hord 1973). 
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Figure 4a shows the temperature profile in the cavity. The strong temperature depression at the leading edge of the 
cavity is evident with the gradual temperature recovery due to condensation in the rear of the cavity. The vapor 
volume fraction is qualitatively compared with a typical flow visualization of the flow (note that the flow conditions 
at which this visualization was done are not reported in the Hord’s report). The overall shape and features of the 
cavity appear to be similar. The computed pressure field indicates strong interaction between the cavity and the 
tunnel wall which is expected due to the relative scales of the geometry. 


Table I. Run Conditions For Liquid Nitrogen Cases. 


RUN NUMBER 

FREESTREAM 
TEMP. (K) 

FREESTREAM 
VEL. (M/S) 

CAVITATION 

NUMBER 

CAVITY 

LENGTH (CM) | 

289C 

88.64 

23.5 

1.55 

2.29 

290C 

83.06 

23.9 

1.70 

1.9 

293A 

77.64 

24.0 

1.75 

1.52 

294F 

77.94 

9.8 

1.78 

1.52 





Vapor Pressure 

Figure 3. Physical Properties Of Liquid Nitrogen (Temperature Dependence). 



Pressure Distribution Cavity Temperature 

Figure 4. Cavitating Hydrofoil Flowfield For conditions of Run 290C in Liquid Nitrogen. 

The quantitative comparisons of pressure and temperature depression in the cavity are compared with 
experimental data in Fig. 5. Note that the pressure values plotted are (P-P vx> ). For a non-cryogenic case this 

value would be zero in the cavity, while values below zero in the cryogenic case indicate pressure depression due to 
thermal effects. In general excellent comparison is obtained for the leading edge temperature depression of 
approximately 2.5 K. The temperature recovery within the cavity compares well both in their slope as well as the 
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length of the cavity. The temperature rise in the cavity closure region shows some differences; the computed 
solution recovers to the freestream value more quickly than does the data, which doesn’t quite fully recover to the 
freestream value. As pa- the discussion by Hord 1 , this probably was due to the unsteady effects in the cavity closure 
region whereby the thermocouples were not always enclosed in vapor giving erroneous readings. The 
instrumentation error given for the chromel-gold thermocouples is 0.20 K and the computed results are within the 
uncertainty of the experiments. 

The comparison of the pressure depression in Fig. 5a also indicates excellent overall comparison with data and 
within the instrumentation error bar of 0.69 N /cm 2 . The leading edge pressure depression is 23.5 percent relative 
to the freestream vapor pressure. It further illustrates why temperature effects have a substantial impact on the 
performance of cryogenic pumps. The experimental data plotted includes both the actual pressure measured 
(symbol: circle) as well as the saturation pressure values (symbol: square) corresponding to the temperature 
measurements. The close match between the actual pressure and saturation pressure values indicates that the 
thermodynamic equilibrium assumption is valid for liquid nitrogen flows. 




Pressure depression. Temperature depression 

Fig. 5. Comparison of Cavitating Hydrofoil Simulations with Data for Liquid Nitrogen. 

V. CRYOGENIC INDUCER SIMULATIONS 

Simulations of a full-scale LOX inducer are presented here and compared with a sub-scale configuration 
designed for water testing. The sub-scale water configuration has a mass flow rate of 250.3 Ibm/s, a rotational speed 
of 4800 rpm, and an inlet water temperature of 85°F. The sub-scale inducer has four blades with an inlet hub 
diameter of 1.654 inches and tip diameter of 6.016 inches. The tip gap clearance is 0.022 inches. The simulations 
were carried out by generating a multi-element unstructured grid for the complete 4-bladed configuration. The bulk 
of the grid cells between the blades is composed of hexahedral cell elements while tetrahedral and prismatic cells are 
used in the tip gap region as well as the region between hexahedral blocks to reduce grid skewness. The total grid 
size was approximately 4.5 million cells and the solution was obtained by decomposing the grid on 32 processors on 
a parallel machine. 

The full-scale LOX inducer was simulated over a range of numbers ranging from 10000 to 30000 to model 
inducer performance from mildly cavitating conditions through head breakdown. Figure 6 shows the vapor volume 
fraction iso-surface ova four different N w numbers. We note the formation of the primary cavitation front at the 
leading edge which grows preferentially near the shroud as the N ss number increases. However, interestingly we 
also evidence the formation of a second cavitation front at a bend in the blade where the blade goes to its full 
thickness. The corresponding pressure distribution on the blade suction side surface is shown in Fig. 7. The loss of 
pressure generation corresponds to the advancement of the vapor formation as expected. The temperature 
depression levels due to thermal effects of cavitation are shown in Fig. 8. Lower temperatures occur in regions 
where vaporization of liquid takes place. This figure indicates that the temperature depression is localized around 
the leading edge of the blade; the primary location of vapor formation, therefore, is at the leading edge while in the 
remainder of the vapor cloud condensation back to liquid is taking place. This is consistent with the attention given 
to the leading edge shape in design tools for inducers. 

The pressure distribution on the both the pressure and suction side of the blade surfaces are shown at different 
Nss numbers in Fig. 9 for both the water and the LOX case at a location close to the blade tip. We observe that the 
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loading profiles for the water and LOX configuration are similar at lower N ss numbers but that the water case breaks 
down earlier. As cavitation sets in, the pressure on the suction side drops to the nominal vapor pressure value and 
this front advances along the suction side consistent with the vapor cloud location. An interesting point to note is 
that while the pressure levels on the pressure side drop under cavitating conditions, the leading edge continues to 
take on load up until an N ss number of 22000. At 25000, the water case has broken down but the LOX configuration 
is continuing to generate head indicating that thermal effects have improved the mean performance of the pump. 

To understand the cavitation loss mechanisms better, we plot the velocity profiles downstream of the blade exit 
for the non-cavitating case as well as the cavitating case at N ss of 25000 in Fig. 10. The axial velocity plots for both 
the non-cavitating and the cavitating case show a slowdown near the shroud due to losses. However, the slow down 
in the axial velocity is more pronounced in the cavitating case indicating larger losses. The differences in the axial 
velocity also lead to differences in the flow angle at the exit. The loss mechanism in the cavitating case was further 
investigated by plotting the streamlines near the shroud (Fig. 11). The streamlines in the cavitating case are affected 
strongly by the blockage due to the vapor cloud. The flow passage gets constricted at the inlet causing the flow to 
accelerate and squeeze through the narrower effective flow passage. Conversely at the cavity closure region, the 
fluid experiences a sudden expansion in flow area as the blockage ends. This causes larger flow turning losses that 
are reflected in the exit velocity profiles. 



Single Phase N ss =20000 N ss =22000 N ss =25000 

LOX Case 

Fig. 6. Vapor Isosurface (0> = 0.4 ) Contours at Various N ss Numbers in LOX Inducer. 



Single Phase ^=20000 ~N SS ~22000 "" ~N SS =25000 

LOX Case 

Fig. 7. Blade Suction Side Pressure at Various Numbers in LOX Inducer. 



Single Phase N ss =20000 N ss =22000 N ss =25000 

LOX Case 

Fig. 8. Blade Suction Side Temperature at Various N ss Numbers in LOX Inducer. 
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Fig. 9. Blade Loading Profiles (Lox Case) at Various Radial Locations in LOX Inducer 


The comparison of the head coefficient generated by the inducer as a function of the N M number are plotted in 
Fig. 12 for both the full-scale LOX case as well as die sub-scale water test. We note that although we did not 
discuss the details of the sub-scale water simulations, the qualitative flow features for the water configuration were 
very similar to the LOX case. In particular, the location and shape of the vapor cloud had similar features although 
the vapor cloud in the LOX case evolved at a slower pace due to thermal depression effects. The effect of the 
thermal depression is apparent on the critical N* number at which breakdown occurs. In the water case breakdown 
occurs in the range of 23000-24000 while in the LOX case breakdown occurs around 30000. Furthermore, in the 
LOX case, the slope of head drop-off is more gradual as the breakdown point is approached probably due to the 
moderating influence of the thermal effects. The single-phase result for the LOX case shows a slighdy lower head 
rise than the single-phase result for water despite the higher Reynolds number for the LOX case. This may due to 
the much higher tip-gap losses in the LOX case; we note that the LOX case was run with the tip gap from the sub- 
scale geometry being scaled up by the same scale factor. 


Axial Velocity 



Tang. Velocity Flow Angle 



Plane 1 

Fig. 10. Velocity Profiles Downstream of Blade Exit Compared for Cavitating and Non-Cavitating Case. 


Single Phase 


Cavitating Case (N ss =25000) 
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Fig. 11. Velocity Streamlines on a Radial Plane Close to the Shroud. 



Nss Number 

Fig. 12. Head Coefficient for Inducer Compared for Water and LOX Case. 

VI. SUMMARY 

A generalized multi-phase formulation has been developed for cavitation in “real” fluid flows (e.g., cryogenic 
fluids) that operate at temperatures close to their critical temperatures. The formulation is based on a thermal 
equilibrium assumption between the vapor and liquid in the cavitating region. Fluid thermodynamic properties are 
specified along the saturation line as a function of temperature using the NIST- 12 data bank. The energy equation 
for the mixture is solved in conjunction with the mass and momentum conservation and the evaporative cooling 
effects of cavitation are accounted for rigorously. 

Detailed validation has been presented to quantify temperature and pressure depression in liquid nitrogen and 
hydrogen by simulating the cavitating hydrofoil experiments of Hord 1 . In general excellent comparison of the 
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leading edge pressure and temperature depression was obtained for all cases. The temperature profiles in the cavity 
closure show more variation and this is in part due to the unsteadiness in the data, which was not characterized. 

The multi-phase formulation was used to simulate a full-scale LOX inducer as well as the corresponding 
subscale configuration designed for water testing. The head breakdown appears to be in the range of 23000-24000 
Nss for the water model while the inducer breaks down at 30000 for the LOX. The increase in suction 
performance is attributed to effect of thermal depression in LOX. Furthermore the head breakdown in the LOX case 
was found to be more gradual than in the water case particularly near the breakdown point The overall cavitation 
characteristics, however, were otherwise similar between the LOX and the water case and in particular a second 
cavitation front was observed at a bend in the blade where it goes to its full thickness. The temperature depression 
was observed primarily at the leading edge with the largest temperature drop being at the leading edge tip. The 
velocity profiles downstream of the blade exit were analyzed. Both the single-phase and cavitating cases showed a 
slow down in the axial velocity near the shroud. However the slow down was much more in the cavitating case due 
to the additional losses in the cavitating case. These additional losses were attributed to the increased flow turning 
losses in the cavitation case due to blockage from the vapor cavity. 
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